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THE STRUCTURE OF GLOBAL ATTRACTORS FOR DISSIPATIVE 
ZAKHAROV SYSTEMS WITH FORCING ON THE TORUS 

M. B. ERDOGAN, J. L. MARZUOLA AND N. TZIRAKIS 

Abstract. The Zakharov system was originally proposed to study the propagation of 



o 

O^ Langmuir waves in an ionized plasma. In this paper, motivated by the work of the first 



and third authors in [H] , we numerically and analytically investigate the dynamics of the 
•^r dissipative Zakharov system on the torus in 1 dimension. We find an interesting family of 

iy~^ stable periodic orbits and fixed points, and explore bifurcations of those points as we take 

^ weaker and weaker dissipation. 

< 

^ 1. Introduction 

a 

i_i In this paper we study the dissipative Zakharov system with forcing: 



(1) 



iut + Uxx + iju = nu + f , x S T = M/(27rZ), tG[0,oo), 

ntt - n^x + 8nt = {\u\'^)xx, 

n(x,0) = uq{x) e H^{T), 

n(a;,0) =no(3;) gL2(T), nt(x, 0) = ni(x) G ^-^(T), f (^ H^{T) 



The original Zakharov system (7 = 5 = / = 0) was proposed in 114J as a model for the 
^ collapse of Langmuir waves in an ionized plasma. The complex valued function u{x, t) 

/\ denotes the slowly varying envelope of the electric field with a prescribed frequency and the 

real valued function n{x, t) denotes the deviation of the ion density from the equilibrium. 
Smooth solutions of the Zakharov system obey the following conservation laws: 

II«(*)IIl2(T) = ||^io||L2(x) 



and 



E{u,n,u){t) = I \dxu\ dx -\ — / n dx -\ — / v dx + I n\u\ dx = E{uo,no,ni) 
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where v is such that nt = Vx and ft = (n + |up)a;. These conservation laws identify 
H^ X L^ X H^^ as the natural energy space for the system. Local and global well-posedness 
in the energy space was established by Bourgain |2|. Lower regularity optimal results were 
obtained by Takaoka in [11 1. The well-posedness theory extends to the dissipative and 
forced system without difficulty [5]. 

In [5], the first and third authors established a smoothing property for the Zakharov 
system, and as a corollary they proved the existence and smoothness of a global attractor 



in the energy space. For a discussion of basic facts about global attractors see 12 and l5|. 
The problem with Dirichlet boundary conditions had been considered in |:6i| and [7] in more 
regular spaces than the energy space. The regularity of the attractor in Gevrey spaces with 
periodic boundary conditions was considered in [JlO]. 

Here, we primarily focus on the dynamics of solutions to (II|. For large dissipation we 
prove that the global attractor is a single point consisting of a unique stable stationary 
solution of the system. Then, we proceed to investigate numerically the case of smaller 
dissipation in the spirit of the numerical exploration of damped-forced Korteweg-de Vries 
equation in |3| and for the Waveguide Array Mode-Locking Model in [13|. In particular, we 
explore equilibrium and periodic solutions, the branching of solutions, bifurcation points, 
period doubling and other interesting dynamical structures that arise. 

The paper is organized as follows. In Section [2| we obtain preliminary estimates on the 
solutions and study the existence and uniqueness of stationary solutions. In Section |3| we 
prove that in the case of large dissipation, the global attractor consists of the unique sta- 
tionary solution. In the remaining sections we study the small dissipation case numerically. 



1.1. Ackno'wledgments. M.B.E. and N.T. were partially supported by NSF grants DMS- 
1201872 and DMS-0901222, respectively. J.L.M. was partially supported by an IBM Junior 
Faculty development award and acknowledges the Universitat Bielefeld as well as the Uni- 
versity of Chicago for graciously hosting him during part of this work. The authors wish to 
thank Jon Wilkening and Roy Goodman for helpful discussions throughout the preparation 
of this work. 



2. Existence of Stationary Solutions and Preliminary Estimates 

We start by obtaining a simple bound on the !?■ norm of the solution. By multiplying 
the ti-equation with u and integrating on T and then taking the imaginary part, we obtain 

d 

This implies by Gronwall's and Cauchy-Schwarz inequalities that for sufficiently large t 
(depending only on the L^ norm of the initial data), we have 

I2 



\u\\l + 2-f\\ug = 2Qjfu. 



(2) ||u||2 < 2 

7 

We now study the stationary solutions of the system (II]). Recall that n is real, and 
throughout the paper we assume that n and rit are mean zero. Let {v, m) be a stationary 
solution of Q. Taking ut = nt = nu = leads to 

v^x + i-yv = mv + f, X e T, 
-m^x = (|f P^ 



J XX: 



The second line of ([3| implies that 771 = — |fp + ax + 6. Therefore, the periodicity and the 
mean zero assumption lead to m = — |vp + 2^||^||2- Substituting to the first equation, it 
suffices to study 



(4) 



9^ ■ 1 „ „2 , ,2 



v = f, xeT. 



Lemma 2.1. Fix f £ L^ and 7 > 0. Any solution f 0/ Ml) satisfies the following a priori 
estimates 

(5) ||Hl2<-||/||2, 

7 

(6) \\vxh < Cmax{j-^f\\lr'\\f\\ll-'^'\\fh). 

Proof. By multiplying Q with v and integrating on T and then taking the imaginary part 
of the equation we obtain that 

7lbll2 = 9 / fvdx. 
This implies ([5| by Cauchy-Schwarz inequality. 
On the other hand, taking the real part we obtain 

1 



Vx\ dx = ||f 11^,4 — — 11^112 — ^ I fvdx. 



By the Gagliardo-Nirenberg and Cauchy-Schwarz inequalities, we have 

||^;x||i<C(||^.||2|b||i + ||f||^ + ||/||2|bl|2). 

This and ^ imply Q. D 

We now prove the existence of an H^ solution f of Q for large 7 and/or small ||/||//i 
which is unique in a fixed ball in H^. Uniqueness in the whole space will follow from 
Theorem 13.11 below. 



Proposition 2.2. Given f G H^{T) if 'y > is sufficiently large, or for given 7 > i/ 
jji is sufficiently small, then we have a unique solution of Q in the ball B := {v : 



Iblbi < ^MhI}. Moreover V G H^{T). 

Proof. First note that by Kato-Rellich theorem the operator ^-^ — 2:^||f||2 + I'Wp is self 
adjoint on L^(T) for v G -L°°(T) C H^{T). Therefore the operator 

^^ 1 „ „2 , ,2 . 

is invertible on L^(T) and we have 

(7) \\R^,v\\l2^l^ < -• 

7 

Let 

T.jiv) := i?-i/. 

It suffices to prove that Tyj has a fixed point in H^. To do that we will prove that T^j is 
a contraction on the ball B. 
By the resolvent identity, 

we have 

Therefore, we obtain 

(8) \\R-!j\\m < Mki + c(7-^/^>7-^/l( - l\\v\\l + \v\')R-!j\\,. 

< ^fc + c(7-^/^)7-^/l -hv\\l + \vn,4R:;!jy 

7 zvr -^ ' -^ 



< 






Note that in the second to last inequahty we used ([T]). Let M = ^||/||hi. The inequahty 
above imphes that for sufficiently large 7 or for sufficiently small ||/||//i, T^j maps B = 
{v E H^ : 1 1 I'll //I < M} into itself. Thus it suffices to prove that T^j is a contraction. Again 
by the resolvent identity and ([8]), we have (for u,v G B) 



. ,|2 II Il2 

\u\\i — f 



ll/'U>|2 l„|2 I NilzMl^ R-1 fll , U-ll2\ 

^^ 1^ + ^ „,i/2 II^IIhij 



71/2 



II, ,9 , ,9 llnllo — llf II9 II llfllH-l 



27r 11^' 72 

Therefore T^j is a contraction on S for small H/Hhi or large 7. Finally by the following 
calculation the fix point v € i? is in H'^iT), 

(9) 11^11^3 = \\T,j{v)\\h-^ < {r')\\f\\m + {r')\\ - ^\\v\\l + MIh^WR^^vIWh^ 

<{r')\\f\\m+C{j-')\\f\\l.r\ 
using the standard elliptic estimate ||(^ +^7) IWh^ < (7 "'^) 11/ II hi- ^ 

3. Attractor in the Case of Large Dissipation 

Recall that the energy space is X = H^ x L? x H^^. We will prove under some 
conditions on 7,5, ||/||//i that all solutions of M converge to the stationary solution 
(v,— |fp + 2^1^112,0) in X as t — )• 00. This also implies the uniqueness of the station- 
ary solution V under these conditions. 

Theorem 3.1. Given ||/||//i and 6 > 0, the following statement holds if ^ is sufficiently 
large. Consider (u(0),n(0), nj(0)) G X where n(0) and nt{0) are mean-zero. Then, the 
solution {u,n,nt) of (II]) converges to the stationary solution {v, — |vp + 2^11^111)0) in X as 
t — )• 00. 



Proof. Given solution {u,n,nt) of Q, let 



{w,z,zt) = {u-v,n+ lup - ^^\\v\\l,nt). 

zvr 

Note that z and zt are mean-zero. The equation for {w, z, zt) is the following 

iwt + Wxx + iiw = z{w + v) — \v\'^w + i^\\v\\\w^ x S T, t£[0,cxD) 



(10) 



|2 hPl 



Ztt - Zxx +^Zt = {\W + V\ 

Fix e > and let 

H = \\dx^{zt + ez)\\l + ||z||^ + 2||u;^.||^ + 2 / z{\w + v\^ - \v\^) + \\w\\l. 

The above quantity H was introduced in |6|| to obtain bounds in the energy space. We note 
that H is bounded by a constant multiple of the energy norm for any fixed e. 
We have 



d , 



= 2 / d-^{zt + ez)d-^[{z + \w + uj^ - \v\'^)xx + (e - (5)zt] 
= -2 /(zt + ez)(z + \w + 7;p - jw^) - 2(5 - e) /"^-^(zj + ez)d-^zt 
= -;^ll-2|l2-2e|k||2-2 / zt(|w + v|^- bl^) - 2e / z(|w + v 
- 2(5 - e)\\d-'ztg - 2e{5 - e) j d-^zd'^zt. 



''-\v?^ 



Using 



we obtain the following energy-type identity 

^ ''d-\z, + ez)||^ + |lklli = -26|kll^ - 2 I z,(|u; + ^1^ - \vf) 

-2e I zi\w + v\^ - \v\^) -{5- e)\\d^^zt 



(11) 



dt' 



1^ l|2 
2 



'l^l|2 



-(5-e)\\d-Hzt + ez)\\'^ + i6-e)e^d-'z 
Now consider the derivative of the remaining terms in the definition of H: 
(12) 2 — ||u;2;||2 = -43f? / w^wt = -'^iW'WxWl - 49 / w^[z{w + v) - \v\'^w] , 



and 



(13) 2— / z{\w + v|2 - Ivp) = 2 / zt{\w + v\'^ - Ivp) + 4K / ztf)^(u' + ?;) 

= 2 / zt{\w + vp - lup) + 49 / zt(}^(t(; + f ) 
— 473? / zw{w + t") + 49 / 2;t;w[|t;| — — - 



l^lli] 



We also have 

(14) 



^tll^lli = ~27||Ti;||2 + 29 / zwv. 



Combining (|11|)-(|14|), we observe 

d_ 
dt' 

+ {6 - e)e^\\d^'^z\\l - 4'y\\wx\\l + 49 / ww^\vf - 47K / zw{w + v) 



-H = -2e\\z\\l - 2e j z{\w + v\'' - \v\^) -{5- e)\\d-^zt\\l -{5- e)\\d-\zt + ez) 



+ 49 / zwit^Nvl 



1 

2^ 



\v\\f^ -27||it;||| + 29 / zwv. 



Hence, 

d 



dt 



H = -eH - e\\zg -{6- e)\\d-'^ztg -{6- 2e)\\d-\zt + ez)\\'^ - (27 - e)\\w 
+ 29 / zwv + {6- e)e^\\d~'^z\\l - (47 - 2e)||u'^||i - 49 / ww^{\v\'^)^ 

— 47K / zw{w + v) + 49 / 2;wtt;[|t;p ll^lli] ■ 

Let e = min(j^, 2;7)- Since z is mean-zero, the choice of e implies 



{6-e)e^d~'zg<'-\\z\\l 



Therefore, we have 
d 



dt 



H<-eH- Uzg - -f\\wg - 2-f\\w^g + 2 



zwv 



+ 4 



WWxi\v\'^)x 



+ 47 



zw{w + v) 



+ 4 



zvw 



27r 



\vg 



< -eH - -\\zg - "f\\wg - 27||u;a;||2 



C 



\z\\2\\w\\2\\v\\m + ll^i'l|2||^«x||2||'y|lH2 + ll-^lbllw^lblbllj/i 



+ 7||z||2||lt'||j:/i||tL' + ^^112 
7 



e^- 7,\\A\l -l\\w\\2-'^l\\wx 



2 



+ [X + XX + XXX + XV] . 
Note that by ^ we have 

7 
for sufficiently large t. Using this, we can bound term XV by 

- lO" "^ e " "^ e " ^112 --^qII 112 -^qII 112 -^qII x\\2, 

provided that e"f » 11/ Hi- Summands X — XXX can be bounded by the same right hand side 
provided that 

^7 > Ibll^i + ll^llffi' ^'^d 7>ll^ll//2- 

By the estimates on v, we see that for fixed 5 and ||/||//i, if 7 is sufficiently large, we have 
for sufficiently large t, 

U < -eH. 

dt ~ 

This implies that H goes to zero as t — )• 00. 
Observe that 



[z{\ 

JT 



W + v\ — \v\ ] 



< C||z||2||'w||i:^i(||'w||2 + ll^^lb), 
and that, by Q and ([s]), (for large 7 and t) we have \\w\\2 + ||w||2 *^ 1- Therefore, we have 

H >C{\\ztfH-. + Ml + \\wfH,). 
This completes the proof. D 

4. Numerical Methods 

In this section we will study the Schrodinger-Dirac model (see, e.g., [5|) numerically. In 
particular, we apply the time-splitting method of |8| as applied to the soliton dynamics in. 



for instance, J9]. Let us recall the equivalent system to ([T]) derived in [s], which are 

{idt + 9^ + i'y)u = nu + f, 



(15) 



{idt - d + i6)n = d{\u\'^), 
{u{x,0),n{x,0)) = {uo,no) e H^ x L? , 



where d = {—dxx)'^- Note, this can be re-written as 



dt 



u 



n 



where 



L 



LI "" I +N{u,nJ), 
n 



idi-j 

-id -6 



N 



and 

-inu — if 
-id{\u\'^) 

The algorithm takes place as a pseudospectral method on the Fourier side, though it im- 
plements integrating factor, time-splitting, fourth-order Runge-Kutta schemes and contour 
integration all at once. The key idea is to look at the evolution over a time step, h, as the 
integral 



Um+l 
nm+1 



^Lh 



+ e 



Lh 



Ur. 



e-'^'N{u{t^ + s), n(t„ + s), /(t^ + s))ds, 



which can be approximated using a Runge-Kutta method (see Cox-Matthews [4]) as 



Urn+l 



^Lh 



+ h-'^L-^ X 



Ur. 



+2 



+ 



4 - L/i + e^'^(4 - 3L/i + {Lhf)] N{um,n^, f{tm))+ 
2 + Lh + e^^{-2 + L/i)] (iV(a^,i, a^,2, /(t„^ + h/2)) 

+N{bm,l,brr,,2j{tn, + h/2))) 

-4 - 3/1 - {Lhf + e^'^(4 - L/i)] iV(c„,i, c„,2, f{tm + h)) 

9 
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7 

Figure 1. A plot of periodic and converged orbits for 7 = .2 to 7 = .45 
with r/ = 1.0 versus the final energy of two distinct branches of converged 
solutions and measured over the average final energy of one period of two 
distinct periodic branches. The solutions were plotted over a choice of ^ = 
-7.0, -4.05, -3.07, -1.0, 0.0, 1.0, 1.95, 3.97, 10.0. 



where 



ttr. 



,Lh/2 



,Lh/2 



Un 



+ L-i(e^^/2 _ Id)N{Um, rim, fitm)), 



n. 



m 



hm = e^"/M "^ \ +L~\e^''/^-Id)N{am,i,am,2j{tm + h/2)), 

Cm = e^'^/^am + L-l(e^'^/2 _ Id){2N{bm,l,bm,2, f{tm + h/2)) - N{Um, rim, f{tm)). 

However, such an algorithm can have problems if L has eigenvalues near 0. To avoid such 
problems the algorithm is slightly modified by evaluating contour integrals over whole discs, 
which are approximated by appropriate Riemann sums. 

5. Numerical Results in the Case of Small Dissipation 

To begin, in an attempt to model the non-trivial dynamics in the Zakharov system, we 

follow some of the ideas in [3] to analyze a series of numerically integrated solutions of 

10 
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Figure 2. A plot of a range of rj from .2 to .8 versus the measured average 
final energy blown-up near a bifurcation point for 7 = .225. The solutions 
were plotted over a choice of A = -4.05, -3.07, -1.0, 0.0, 1.0, 4.0. 
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Equilibrium Branch 1 
Periodic Branch 2 




Figure 3. Plot of the Energy versus (appropriately translated) time for 
initial A = —4.05 (Nearly Doubly Periodic), —4.0 (Periodic Branch 1), -5.0 
(Periodic Branch 2), 0.0 (Equilibrium Branch 1) for 7 = .25 and r/ = 1.0. 
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Figure 4. Plots of the Dirac Final Energy versus Schrodinger Final Energy 
for two solutions that converge to an equilibrium point of the dynamics 
with 7 = .25, r] = 1.0 (left, Equilibrium Branch 1) and 7 = .35, rj = 1.0 
(right, Equilibrium Branches 1, 2). Specifically, we take A = 0.0, —1.0 and 
A = —4.0, 0.0 respectively and plot the parametric evolution from t = 20.0 
to i = 40.0. Note that the trajectories converge to the same equilibrium 
point for 7 = .25 and to different equilibrium points for 7 = .35. 
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Figure 5. Plots of the Dirac Final Energy versus Schrodinger Final Energy 
in the case of periodic orbits in the dynamics with 7 = .25, r] = 1.0 and A = 
—5.0 (left. Periodic Branch 2), A = —4.05 (right. Near Doubly Periodic). 
We plot the parametric evolution from t = 20.0 to t = 40.0. The solution on 
the right will, along upon further evolution, eventually converge to Periodic 
Branch 1 as in Figure |3j The doubly periodic dynamics occur only for a very 
small window of A values, A G [—4.05, —4]. 
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Figure 6. Plot of the Dirac Final Energy versus Schrodinger Final Energy 
in the case of a stable doubly periodic orbit from Periodic Branch 2 corre- 
sponding to 7 = .225, r/ = 1.0, and A = 1.9005. We plot over the parametric 
evolution over two full orbits from t = 180.0 to t = 200.0. The doubly 
periodic dynamics occur only for a very small window of A values, A £ 
[1.9,2]. 



(15). In our numerical experiments we observe a great deal of energy exchange between 



the Schrodinger and Dirac solutions, hence we will focus on relatively small energy initial 
data in order to justify that our numerics are valid on long time scales. If the Fourier 
modes become too large at the edges of the spectrum, we do not consider the solution to 
be appropriately accurate, hence all simulations included here will have small contributions 
at high frequency. The time scale on which we integrate is generally T = 50.0 with the 
time step /i ~ le — 4. We will begin by taking 32 Fourier modes on which to evolve. In 
addition, our contour integrals in the numerical evaluation of L^^ will be taken as a mean 
of 64 equidistributed points along the disc. 
We look for solutions with 



u{0,x) = Asm{x), n(0,x) = 0, (5 = 7, / = ??sin(x), 

13 




Figure 7. Plots of the Energy versus (appropriately translated) time for 
periodic doubling along Periodic Branch 2 of solutions with 7 moving from 
.25 (Period 1) to .225 (Period 2) with A = -5.0 and ^ = 1.90 respectively. 



where we vary 
(16) 



10<^<10, .2<7<1.0, .2<r?<1.0. 



For the range of rj in (16), given 7 > 1, we observe that the dynamics tend to a fixed 



equilibrium solution for a large set of A values, putting us solidly in the dynamics of Section 
[3} However, for 7 much smaller, we observe much richer dynamics in the phase space, 
particularly in the form of quasiperiodic orbits and equilibrium solutions. We use the phrase 
quasiperiodic orbits to mean that the orbits are periodic up to preset numerical integration 
parameters. Namely, we assume that the solution returns to a previous configuration within 
an error of at most 5e — 3. However, many such orbits can then be integrated over many 
periods and are in that sense quite numerically stable. 

As one may expect, the nature of our orbits can change as we vary 7 or r/, which accounts 
for movement in the bifurcation diagram presented in Figures [T] and [2] The other figures 
present particular solutions for rj = 1 and for various values of 7 and A justifying the 
bifurcation diagram in Figure [T| In Figure [3j we present three solutions corresponding to 
periodic branches and the equilibrium branch from Figure [T] for 7 = .25, and an unstable 

14 



near doubly periodic solution that eventually converges to the solution in periodic branch 
1 (also see Figures Is] and p| . In Figure El we present sample solutions corresponding to 
equilibrium branches. For 7 = 0.25, we have one equilibrium solution and for 7 = 0.35, 
there are two distinct equilibrium solutions. Finally in Figure [7j we present period doubling 
by plotting the time behavior of two limiting solutions from Periodic Branch 2 with 7 moving 
from 0.25 to 0.225. 

6. Discussion 

We have analytically and numerically observed rich dynamics in the dissipative periodic 
Zakharov system with forcing. Open problems for future consideration include understand- 
ing the large exchange of energy from Schrodinger to Dirac, classifying dynamics for a larger 
range of energies, finding more bifurcation points, etc. For small 7 values, we observe nu- 
merically that there is a great deal of energy transfer from the Schrodinger equation into 
the Dirac equation at the outset of the dynamics. As a result, it is a challenge to probe 
extremely large energies within the current numerical framework. In addition, such energy 
transfers seem to dissipate through the damping on incredibly slow scales, making it diffi- 
cult to search for meaningful dynamics such as further period doubling, Hopf bifurcation 
branches and more as the orbits behave do not settle onto a fixed point or a stable orbit 
for very long time scales. However, we mention that using the orbits we have shadowed 
here, it is likely that a variation of the gradient methods in [1] could allow one to construct 
nearby periodic solutions with great accuracy and hence move along the solution branches 
in a more robust manner. 
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